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Abstract: Traditional object tracking technology usually regards the target as a point 
source object. However, this approximation is no longer appropriate for tracking extended 
objects such as large targets and closely spaced group objects. Bayesian extended object 
tracking (EOT) using a random symmetrical positive definite (SPD) matrix is a very 
effective method to jointly estimate the kinematic state and physical extension of the target. 
The key issue in the application of this random matrix-based EOT approach is to model the 
physical extension and measurement noise accurately. Model parameter adaptive 
approaches for both extension dynamic and measurement noise are proposed in this study 
based on the properties of the SPD matrix to improve the performance of extension 
estimation. An interacting multi-model algorithm based on model parameter adaptive filter 
using random matrix is also presented. Simulation results demonstrate the effectiveness of 
the proposed adaptive approaches and multi-model algorithm. The estimation performance 
of physical extension is better than the other algorithms, especially when the target 
maneuvers. The kinematic state estimation error is lower than the others as well. 

Keywords: extended object tracking; Bayesian approach; random matrix; interacting 
multi-model algorithm; model parameter adaption 
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1. Introduction 

During the past several decades, the resolution of radar and other sensors continues to increase, 
thereby generating complex requirements for object tracking (OT) technology [1]. Traditional OT 
technology characterizes the target as a point source object whose physical extension is ignored [2]. 
The centroid is usually utilized to represent the target, and its kinematic state such as position, velocity, 
and acceleration is estimated in OT. This simplification is meaningful only when physical extension 
can be neglected compared with the sensors' measurement error. However, considering the increasing 
sensor resolution, physical extension state should be estimated jointly in OT approaches, especially 
when tracking relatively large targets such as aircraft carriers and large transport planes. This method 
is referred to as extended object tracking (EOT) [1,3-5]. When targets form a closely spaced formation 
(e.g., an aircraft formation), the individual target becomes indistinguishable from the formation 
because of limited sensor capabilities and sensor-to-target geometry. Thus, closely spaced group 
objects can also be considered extended objects [1,3]- Data association in EOT is difficult because 
extended objects produce a highly varying number of measurements corresponding to a single extended 
object [3]. 

A number of techniques and approaches have been proposed to solve the EOT problem, which 
requires estimating physical extension and the centroid' s kinematic state jointly. Baum et al. [6] 
proposed a novel modeling approach to describe physical extension. The proposed model was called 
random hypersurface model. A recursive Bayesian estimator was derived based on this model. 
Mahler [7,8] apphed probability hypothesis density (PHD) filter and cardinahzed PHD (CPHD) filter 
to EOT, including both extended and unresolved targets. Vo et al. [9] presented a mathematically strict 
Bayesian filter for EOT, which could be reduced to a CPHD filter under some assumptions when 
tracking a single target. Carmi et al. [10] utilized the Gaussian mixture model to describe the 
time- varying dynamics of group objects. A Markov chain Monte Carlo particle filter for multi-target 
tracking was then obtained. Multiple hypotheses tracking was also a feasible option for EOT [11]. 
Zhang and Bar-Shalom [12] proposed a new framework based on multiple kernel centers for visual 
object tracking. 

Koch [1] proposed the random matrix-based approach for extended object tracking. A recursive 
form filter is derived within the Bayesian framework, and the final form of this random matrix-based 
filter (RMF) is more or less similar to that of the standard Kalman filter. RMF describes the physical 
extension of the target with an ellipse or ellipsoid. Thus, RMF is reasonable and can be applied in 
many practical applications as shown in Figure 1. Given that symmetrical positive definite (SPD) 
matrices can represent all the ellipses/ellipsoids centering at the origin, RMF utilizes the random SPD 
matrix to characterize the physical extension of the target. The random SPD matrix is assumed to obey 
Wishart-related distributions. Then the Bayesian formula can be applied to obtain a very simple filter 
according to the properties of the Kronecker product and Gaussian and Wishart-related distributions. 
Feldmann [3] introduced a more accurate distribution of measurement noise with respect to the target's 
physical extension and derived a new RMF approach. This new approach can be integrated into an 
interacting multi-model (IMM) framework to improve estimation performance if the object state 
switches between maneuvering and non-maneuvering. However, the derivation process cannot be 
maintained within the Bayesian framework, and several approximations are applied to obtain a 
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recursive form filter. Lan and Li [4,5] further improved extension dynamic model and measurement 
noise model by describing extension evolution and observation distortion with real matrices Ak and Bk. 
This improvement is logical given the relationship between the SPD matrix and ellipse/ellipsoid. It can 
be applied within the Bayesian framework, which requires only an insignificant approximation. 

Figure 1. Description by ellipses of the target's physical extension. A large ship (left) and 
an aircraft formation (right). 





RMF approach is much simpler than the other EOT approaches mentioned above, thereby making it 
more promising. The simulation results demonstrate the effectiveness of the three abovementioned 
RMF approaches. However, the physical extension estimation error obtained by these approaches 
significantly increases when the target maneuvers (e.g., turning motion). The extension dynamic model 
in [4] is improved in this study based on RMF approach to simplify the description of extension 
evolution. Model parameter adaptive approaches for both extension dynamic and measurement noise 
are proposed based on properties of SPD matrix and ellipse/ellipsoid. The proposed adaptive RMF is 
called ARMF. The preliminary results of single-model ARMF have been published in a previous 
conference paper [13]. A multi-model algorithm is proposed by integrating ARMF into the widely 
utilized IMM framework to further improve the performance of RMF when the target maneuvers. 
Simulation results show that the ARMF approach and ARMF-based multi-model (ARMF-MM) 
algorithm are effective. 

2. Bayesian EOT Using Random Matrix 

The dynamic state of an extended object at scan k is described by both the centroid's kinematic state 
Xk and physical extension Xk, where Xk is an ^'J x 1 random vector and Xk is a d x d SPD random 
matrix. Spatial dimension d is usually set to 2 or 3 in EOT. s is the dimension of the state in one spatial 
dimension, s = 2 means that position and velocity are considered. If acceleration is also considered, 
then s = 3. 

The basic idea of random matrix-based EOT is to characterize joint density p[xk,Xk\Z^] as a product 
of a Gaussian density and an inverted Wishart density based on Bayes' formula: 

p[Z,\x„X,]xp[x„X,\Z'-'] 



p[x„X,\Z'] = 



llp[Z,\x„X,]xp[x„X,\Z'-']dx,dX, 
= -p[Z,\x„X,]xp[x„X,\Z'-'] 

C 

where - i^lVrU denotes the set of rik measurements at scan nk> I is assumed to be independent 
of Xk and Xk, ^ - {ZjYj=i , and c always denotes the normalization factor in this paper. 
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This section summarizes the key results and steps. More detailed assumptions and derivations can 
be found in [1,3,4,14]. 

2.1. Dynamic and Measurement Models 

The dynamic model of Xk for EOT is very similar to the models utilized in many Kalman-related filters: 

Xk=Pk^k-i^^k (2) 
where Fj^=Fj^0I^, "(8)" is the right Kronecker product [15], represents the state matrix in one 

spatial dimension, and Id is the J x J identity matrix. Independent process noise follows a normal 
distribution [1]: 

co,-mA®I,) = -^(0,Q,) (3) 

where N(/i,2) stands for normal distribution with mean /# and covariance matrix i7, and Dj^ denotes the 
covariance matrix of cOk in one spatial dimension. Fj^ and Dj^ can be selected according to specific 

conditions of EOT. Singer's model [1], constant velocity model, and constant acceleration model [16] 
are all feasible alternatives. 

The following assumptions are introduced for physical extension [1,4]: 

p[X, I Z,_J = W(X,;S„A,X,_,A^ Id,) (4) 
p[X,_,\Z'-']^TW{X,_-v,_,,_,,X,_,,_,) (5) 

where dk> d - \ stands for the degrees of freedom, Vk-i\k-i > 2(i is a scalar parameter, subscript 
^-1|^-1 indicates the estimated value of corresponding variable at scan A: - 1, and A is a J x J 
nonsingular real matrix describing physical extension evolution. 
W(Z; a,A), the Wishart density of SPD matrix Z, is defined by: 

^A)^v A\ 1 I . |-(l/2)« I ^|(l/2)(a-^-l) 1 ^-Ivn 

W{X\a,A) = -\A\ \X\ etr[ — A Z] (6) 

c 2 

where etr[-] is short for exp[trace(-)], a>d, and expectation E[X| = aA. 
IW(Z; a,A), the inverted Wishart density of SPD matrix Z, is defined by: 

XW(X',a,A) = -\A\ \X\ etr[ — AX ] (7) 

c 2 

with a>2d + 2 and expectation E[Z] = A/ (a -2d -2). 

The assumption in Equation (4) is different from that in Equation (10) in [4]. It is simpler and more 
plausible, which will be explained in Section 3.1. 

Measurement is modeled as follows: 

zl=H,x,+v[ (8) 

where H^^ =Hj^<S)I^ and Hj^ = [1,0, • • - ,0] denotes the meas-urement matrix in one spatial dimension [1]. 

Gaussian measurement noise v[ is independent of vi(/ = \,2,...,nk,j ^r). Its distribution is assumed 
to be [4]: 
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v[~M{0,B,X,BD (9) 

where Bk is a d x d nonsingular real matrix describing the measurement distortion of extension. 
Equation (9) indicates that the measurements are affected by physical extension Xk. 

2.2. Prediction 

The prediction density of Xk and Xk can be factorized as: 

p[x,,X, I Z'-'] = p[x, I X,,Z'-']xp[X, I Z'-'] (10) 

The first factor on the right of the previous equation is assumed to have the following structure [1]: 
p[x, I X,,Z'-'] = j p[x, I X,,x,_,,Z'-']xp[x,_^ I X,,Z'-']dx,_, 

= lmx,;F,x,_,,Q,)xX(x,_, )d^.-. (11) 

where subscript k\k-l indicates the predicted value of corresponding variable at scan k . 
We then obtain the prediction equations of kinematic state: 



^k\k-i ~ ^k^k-i\k-i 



(12) 



P =F P +0 

p[Xic\Z'^~^] can be calculated based on Equations (4) and (5): 

p[X, I Z*-i] = J p[X, I X,_Jx p[X,_, I Z'-']dX,_, 

= j WiX, ; S, , A,X,_X IW(X,_, ; v,_,,_, , X,_,,_, )dX,_, ( 1 3) 

= QB^ iX,;SJ2, (v,_.,_. -d-V)/2,A, X,_,,_X 1 5„0) 

where GB "(O is the probability density function (pdf) of "generalized beta type IF' distribution. An 
inverted Wishart distribution is utilized to obtain a recursive-form estimator and approximate GB 
via first and second moment matching [1,4]: 

p[X, I Z^-'] « I>V(X,;v,,_i,X,,_,) (14) 

where: 

_ ^,(V,H-l)(2,,.-l)(V.-2) 

and: 

^k\k-i = (^kik-i -^d- 2)A,X,_i|,_iAj / (16) 

with Xk-i = Vk-i\k-i -2d -2. 

The substitution of Equations (11) and (14) into Equation (10) yields: 

p[x„X, I Z'-'] = A/'(x,; j:,„_i,P,„_,)x2>V(Z,;v,„_„Z,„_,) (17) 
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2.3. Update 

In consideration of Equation (1) and based on Equations (8) and (9), likelihood function p[Zk\xk,Xk] 
is obtained as follows [4]: 



p[Z,\x„X,]ccp[Z,\n„x„X,] = Ylm:;H,x„B,XX) 

r=l 

X Af(z, ,H,x„ ^J^^)W(Z, ; - 1, B,X,Bj) 



(18) 



where: 



Zk= — llz[. Z^ = Y^{z[-Zk){z[-Zkf (19) 



The substitution of Equations (17) and (18) into Equation (1) yields [1,4]: 

p[x„X, I Z^] X M{z,;H,x„ ^'^^'^h mZ,;n, -l,B,X,B^) 

(20) 

)xIW(X,;v,,_i,Z,,_i) 

The product of the two Gaussian probability density functions in the previous equation can be 
converted to [1,4]: 

B X 

^i^k ; -^.1.-1 ' X ^i^k '^H,x,, ' ' ' ) = Af(z, ; H,x,,,_, , ) x Af(x, ; x„„ J (2i) 



"k 



with: 



and: 



^k\k ^k\k-i ^k (^k ^k^k\k-\) 
^k\k ~ ^k\k-i ~ ^k^k\k-i^k 



C -TIP HJT , ^k^k^k ^ ZI U njT . ^k^k\k-\^k 



(22) 



k\k-l '■^ k'- k\k-\'-^ k ' " A:' A: 

^k (23) 

^ Pkik-i^k ^k\k-i 

where Xk\k-i is used to approximate Xk. 

Based on Equations (21) and (20) can be rewritten as [1,4]: 

p[x„X,\Z']ccW(Z,;n,-l,B,XX)^(X,;v,,_„X,,,_,) 



(24) 



Subsequently, the first three pdfs on the right side of the previous equation can be combined as 
TWiXk, Vk\k,Xk\k) with: 
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(25) 



where X]^l_^ and Sj^ll\ are square roots of SPD matrices Xk\k-\ and S^.^.^ , which can be obtained by 

Cholesky factorization. They are utihzed to keep Nk\k-\ and Xk\k in positive definite or positive 
semi-definite structure [3]. The derivation of Equation (25) is similar to the corresponding process in [1]. 

Then, p{xk,Xk\Z^^ is converted to the product of a Gaussian density and an inverted Wishart density 
as follows: 



where ak\k - yk\k - 2d-2> 0. 

The Bayesian recursive estimator for EOT with random matrix (i.e., RMF) is composed of 
Equations (12), (15), (16), (22), (23), (25) and (27). A notable similarity exists between the standard 
Kalman filter and the centroid's kinematic part of RMF. 

3. Model Parameter Adaptive Approaches 

For a certain model used in multi-model (MM) algorithm, the model parameters include in 
Equation (3), dk and Ak in Equation (4), Bk in Equation (9), etc. The model parameter adaptive 
approaches act on extension dynamic model parameter Ay^ and measurement noise model parameter Bk 
based on the relationship between the SPD matrix and ellipse/ellipsoid. The basic principle of 
parameter adaptive approach when d = 2 and d = 3 is the same, d = 2 is regarded as an example to 
introduce the parameter adaptive approaches. 

The following assignments are adopted in the sequel for the purpose of convenience: d = 2, s = 3, 
and Hj^ =[1,0,0], which means that the kinematic state in one spatial dimension is [position, velocity, 
acceleration]^ and only the position coordinates of the target's centroid are measured. Therefore, the 
ellipse describing physical extension can be represented by [4]: 



with X denoting the coordinates of the points on the ellipse. 

The following factors must be identified when determining an ellipse on the 2D plane. 

(1) Location: the coordinates of the ellipse's center point, which are determined by 
kinematic state Xk, 

(2) Size: the lengths of the semi-axes, which are equal to the positive square roots of XkS 
eigenvalues; 



p[x, ,X,\Z']cc Mix, ; P,, )TW{X, ; X,, ) 



(26) 




(27) 



{x-H,x,fXl\x-H,x,) = \ 



(28) 
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(3) Orientation: the direction of the semi-axes, which is determined by the angle between 
either of the semi-axes and either of the coordinate axes when d = 2. Several 
definitions of orientation angles are available. For example, the angle between the 
major semi-axis and x axis is suitable when the counterclockwise direction is 
considered the positive direction. 

Size and orientation are determined entirely by physical extension Xk. As long as the size or 
orientation of an ellipse changes, the SPD matrix representing the ellipse would change. This 
occurrence means that the SPD matrices represent only the ellipses centering at the origin. However, 
all the ellipses centering at any position on the 2D plane will be covered if Xk is combined with 
kinematic state Xk. 

The following lemma is necessary to decompose SPD matrix Xk into a form wherein the 
relationship between Xk and the ellipse can be easily analyzed. 

Lemma 1. If X is a d x d SPD matrix, then X can be transformed to a diagonal matrix by up to 
— ^ ''Givens rotations [1 7]. When d = 2, X can be transformed as follows: 



where R((p^) = 



smg? 

G 



G 



sm (p cos (p 



R\(p'')XR((p'') = dmg(a;,a;) (29) 
is the rotation matrix, and g\ are the eigenvalues of X , 



(D G [ , — ) , and 

2 2 



qP = arctan [^p± <^1 + j 



(30) 



with p = 



Z(l,l)-Z(2,2) 



2Z(1,2) 



Therefore, 



Z=/f(^^>°)diag(a-;,a-;)/f^(^2>°) 



(31) 



Proof of Lemma 1. Expand the left portion of Equation (29), then the following equation is 
obtained by setting the non-diagonal elements to 0: 

Z(l, 2) tan^ (<^° ) + (Z(l, 1) - Z(2, 2)) tan(^° ) - Z(l, 2) = 0 (32) 

Z(l,l)-Z(2,2) 



Let P denote 



2Z(1,2) 



; thus: 

tan(^°(Z)) = -p±^^^ 



(33) 



Q.E.D. 

Decomposing SPD matrix X into the product of orthogonal matrices and a diagonal matrix is 
relatively easy in mathematics. Except for "Givens rotations," singular value and eigenvalue 
decomposition are also possible. However, the orthogonal matrix contains not only rotation 
transformation but also symmetry transformation (rotation matrix is the orthogonal matrix whose 
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determinant is +1). The form and sequence of symmetry transformation are uncertain. Thus, neither of 
the two decomposition methods is apphcable for model parameter adaption. 

According to Lemma 1, SPD matrix Xk representing a 2D eUipse could be factorized as follows: 

a; = B(^p^)diag(a^i.a'^,)J?^(^z?,) 

^ R(<p,)V,R\9.) ^^^^ 

where cr^j and cr^^ eigenvalues describing the ellipse's size, and e 



n n 
■2'2 



is the 



orientation angle describing the ellipse's orientation. 

3.1. Extension Dynamic Model Parameter Adaption 

From the extension dynamic model in Equation (4), we obtain: 

E[XJX,_J = A,X,_A' (35) 

No dk exists on the right side of Equation (35), which is different from Equation (11) in [4]. 
However, employing Equation (35) allows A], to cover all the possible changes in the size and 
orientation of the ellipses represented by SPD matrices X]^ and X]^.\. According to Equation (15), 8^ is 
still incorporated in the Bayesian RMF estimator, describing the uncertainty of physical extension 
evolution [1]. Therefore, the improvement of Equation (4) is reasonable and meaningful. 

For parameter adaption when J = 2, we propose to calculate Ay^ by the following equation: 

A,^R{etw^g{rt\rt') (36) 



where 0. e 



\ 

describes the change in ellipse's orientation and /^'^>/^'^>0 are scalars 



TT n 

describing the change in size. 

Although the range of A is reduced, Ak defined by Equation (36) is still able to approximately 
describe all changes in the ellipses' size and orientation. 

According to Equations (35) and (36), 6^ describes the change amount of orientation angles of the 
two ellipses represented by Xk and Xk-\ respectively. Then 0^ could be set to Is^cpk = (pk~ (Pk-u where 
(pk and (pk-i are orientation angles described in Equation (34). Let (p^{Xk) and (p^{Xk-\) denote the 
angles calculated by Lemma 1. If the difference between ^^(•) and is a constant 
{i.e., (p^iXk) - (p^{Xk-\) = Ag)k), then 6^ could be set to (p^iXk) - (p^{Xk-\). However, no matter 
or in Equation (30) is selected, the difference between ^^(•) and cpk is not a constant, i.e., 
(p^(Xk) - (p^{Xk-\) ^ /^cpk' Thus, the following two restrictions are introduced. 

First of all, at the very beginning of EOT, the initial value of extension state should be selected by: 

^0 =^(%)diag(a-oi,cr,2)/f^(^o) (37) 

According to the previous discussion and Equation (34), X^ should be able to cover all the possible 

SPD matrices representing the ellipses centering at the origin. 

Secondly, we select or ' -" in Equation (30) to satisfy the following condition: 

sign(cr^ - cr^) = sign(cro ^ - cr^ 2) (38) 
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where cji and 02 are the same as those in Equation (29), and sign(-) is the sign function. 

If these two restrictions are met, then ^^(^ differs from cpk by a constant during the whole EOT 
process, i.e., qP{Xk) - (p^(Xk-i) = Acpk). This property is verified by the simulation results. 

However, Xk and Xk-i are unknown and Xj^ is also unavailable at the beginning of period k, thus 
we use Xj^_^ and ^^_2 to approximate Xk and Xk-i. Ot is then calculated according to the change of 
orientation angles of the ellipse represented by ^ and X^_^ , 

From the above conditions, 6^ can be calculated by: 



eigenvalues of X^_^ and X^_^ . When the size of the object does not change rapidly and sharply, the 



where tr(-) is short for trace(-). 

Therefore, Ak can be adaptively determined according to Equations (36), (39) and (40). 

3.2. Measurement Noise Model Parameter Adaption 

Let 21 denote the real co variance matrix of measurement noise v\. According to Equation (9), 
BkXkB\ is used to approximate iT^. Thus, the basic idea of Bk adaption is adjusting Bk to meet 
BkXkB\ = 2^ as accurately as possible. Given that Xk and iT^ are unknown and is also unavailable at 
the beginning of period k, Bk is calculated according to the difference between ^ (k > I) and X. 
where X is an approximation of i7^. Let Wk denote the noise covariance matrix of measurement sensor, 
which is usually known as prior information. Then X is calculated by [3]: 




4 * 

following /'i^ ' calculation method is suggested. 




(40) 



(41) 



X 



k\k-i 



where // is a scaling coefficient, X^|^_jis utilized to approximate Xk, and X^ 



-2d -2) 



klk-l 



according to Equation (14). 

Similar to Equation (36), Bk is calculated by: 






(44) 
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33, Model Parameter Adaption When d = 3 

For 3D EOT (d = 3), Xk can be transformed to a diagonal matrix by up to three "Givens rotations", 
which means that at most three angles are required to determine a 3 x 3 SPD matrix [17]. Thus, the 
following equation is introduced to calculate Ay^: 



be noted that the three rotation matrices in Equation (45) must be maintained in the same order during 
the whole EOT process. In addition, when utilizing Lemma 1 to decompose Xk, the order of rotation 
matrices should be the same as that of Equation (45). The adaptive approach of Bk when J = 3 is 
similar to that of Ak. 

4. Multi-Model Algorithm 

The proposed ARMF approach is integrated into a multi-model algorithm framework [18] (e.g., the 
widely used interacting multi-model (IMM) algorithm [19]) to further improve the estimation performance. 

4.1. Moment Matching 

Moment matching is applied in step initialization and fusion of IMM. Unlike the standard IMM 
algorithm [18], moment matching in ARMF-MM algorithm involves both kinematic state Xk and 
physical extension Xk, this condition makes the moment matching method more complex. 

Step fusion is regarded as an example in this section to introduce moment matching method. In the 
first place, the substitution of Equation (27) into Equation (26) yields: 




p[x, ,X,\Z']cc Mix, ; X,, , P,, )2W(X, ; v,, , a,^,X, ) 



(46) 



Moment matching in step fusion requires that the following approximation be solved: 



k ' ^k\k ' ^k\k ^■^ k\k / 

k ' ^k\k ' ^k\k 



(47) 



where j = 1,2. ...,n denotes model j and n is the number of models. 

The moment matching method for kinematic state Xk is similar to standard IMM: 

-^kik j^i^k-^kik 

where /-4. - pb^d 1 and indicates that model j is in effect at scan k. 



(48) 
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The problem of moment matching method for physical extension is that the method considers Xk\k 
and Vk\k in Equation (25) simultaneously. The mean square error e^. of extension estimation is first 
defined as follows [3]: 

e, ^ tr(E[(.Y, - X,f \Z']) = ti-(E[X; | Z*] - A';) 



K*+i)K*-2) 

-2) 

Expand Equation (49), and the following equation is derived [3,20]: 



(49) 



+ 2d + 2 (50) 



where: 

— \\2 



ft=4+(tr(Z,)) +i,(xl) 
^,=2[^,+tr(X,^)] 



(51) 



Detailed derivation and explanation can be found in [3,20]. Then the moment matching method for 
physical extension is as follows: 

(52) 

^.=EM/'/(^/+tr((x/-x,)^)) 

After and are determined, Vk\k is calculated based on Equations (50) and (51). Thus, Xk\k is 
obtained as follows: 

= ^m^k = (^kik - 2^ - 2)^, (53) 
4.2. ARMF-Based Multi-Model Algorithm 

The four steps of ARMF-IMM algorithm are summarized below: 
(1) Reinitialization 

Let n'^ denote the probability of transiting from model / to model j. n\\ is the probability that 
model / is in effect at scan k-\ given ml and ^ 

^ p[mi\m[_„Z'-'}pH_,\Z'-'] 

p[ml I Z'-'] (54) 

i, i i 

c 

Reinitialization of each model is achieved by replacing corresponding probability variables in 
Equations (48)-(53): 
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pO] ^'V (p'- +/^' )(.)'] 

^l\=T" <vy; 1=7" I I 

t-i Lui-x k-i k-i /:.w(^,,._^^_^_2j_2j (55) 

(2) Model Parameter Adaption and Filtering 

For model j, j = 1,2....,^, calculate A^^ according to Equations (36), (39) and (40). Then the RMF 
approach as summarized in Section 2 is performed and Bj, is calculated according to Equations (42)-(44) 
for model j. 

(3) Probability Update 

According to Equation (54), probability /4 is calculated by: 

/4 =pH I z']=-p[z, I ^<,z'-']pH I z"''] 

C 

^-^Hp[mi\Z''-'] (56) 



c 

Likelihood A^^ is derived as: 



= II KZ, I mi , X, , XMx,.X, I mi , Z'-']6x,dX, 

where p[xk,Xk\m^j^,Z^~^] is calculated by Equation (17) and p[Z^\m^j^,Xk,Xk] by Equation (18). 
The substitution of Equations (17) and (18) into Equation (57) approximately yields: 

r.[K+v4_,-d-2)/2] 



(57) 



r.[(v4_,-rf-i)/2] 



(58) 



xBiXi,_,(Biy 2 x5/X4_,(5/)'+Z, 



where ^ 4_i is used to approximate X{ in S/i^,^ , and Td['] is the multivariate gamma function. The 

derivation of A is quite similar to that of Equations (11) and (13). When rik = 1, 
Ai=M(z,;H,xi,,_„Si,_,). 

(4) Fusion 

The moment matching formulas of fusion are provided by Equations (48)-(53). The final result of 
extension estimation is defined in Equation (27). 
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5. Simulation Results and Discussions 

Two typical EOT scenarios are simulated in this section to compare the proposed ARMF approach 
and other RMF approaches. 

5.1. Scenario 1 

The route of an aircraft carrier, which is approximately 300 m long and 80 m wide, is shown in 
Figure 2. The route starts from origin (0, 0)^, moves along the trajectory, and then maintains a uniform 
circular motion on the 2D plane. Constant velocity is set to 27 knots (approximately 50 km/h). The 
aircraft carrier's extension is characterized by an ellipse with radiuses 170 m and 40 m in RMF 
approach. Scan period T is set to 10 s. The measurement scattering centers are assumed to be 
distributed uniformly over the physical extension {i.e., the ellipse) [3], and the number of 
measurements in every scan obeys Poisson distribution with mean 20. The sensor measurement noise 
follows a zero-mean Gaussian distribution with covariance matrix = diag([50^,20^]) 

Figure 2. Target trajectory of simulation scenario 1. 
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The single-model ARMF is compared with Koch's approach [1] and Feldmann's approach [3] 
through A/m = 500 Monte Carlo runs. 

(1) Koch 's approach', the temporal decay constant x is set to 87; 

(2) Feldmann 's approach: Q.25Xk + Wk is utilized as the measurement noise covariance matrix; and 

(3) ARMF: rj in Equation (42) is set to 0.25 via moment matching [3]. 

Average 0^ , which describes the change in the ellipse's orientation angle, is shown in Figure 3. 
The curve shows that the model parameter adaptive approach is effective and that ARMF is able to 
adjust according to the extension's changes. 

The root mean square error (RMSE) of physical extension Xk is defined by [3]: 



(59) 



where stands for the expectation of physical extension in Ith Monte Carlo run. 
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Figure 3. Average 0^ in simulation scenario 1. 
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The RMSE of physical extension of the three RMF approaches is shown in Figure 4. When the 
object makes a turn (uniform circular motion after scan 10), the extension RMSE of ARMF becomes 
lower than that of the other two approaches. 

Figure 4. RMSE of physical extension in simulation scenario 1. 
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5.2. Scenario 2 

The extended object scenario in [3] is simulated below. The aircraft carrier whose size is the same 
as the one in scenario 1 starts from the origin and moves along the trajectory on the 2D plane as shown 
in Figure 5. Three turns are implemented. 

The following three approaches are compared through A^m = 500 Monte Carlo runs. 

(1) Koch's approach: parameters are the same as those in scenario 1; 
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(2) Feldmann 's MM approach: parameters are designed similar to those in [3]; and 

(3) ARMF-MM: three models are used. Model 1 with low kinematic process noise and low 
extension agility (small (g^ , large Sj^ ), model 2 with high kinematic process noise and high 

extension agility (large (g^ , small Sj^ ), and model 3 with moderate kinematic process noise and 

high extension agility (moderate , small Sj^ ) are utilized, rj = 0.25 for all three models. 



Figure 5. Target trajectory of simulation scenario 2. 
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The simulation results reconfirmed the validity of the model parameter adaptive approach. ARMF-MM 
identified the angle changes correctly when the object maneuvers as shown in Figure 6. 

Figure 6. Average of model 3 in simulation scenario 2. 
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The position RMSEs of three RMFs are close to one another. The error curve is shown in Figure 7. 
The position RMSE of ARMF-MM is lower than that of the others. The RMSEs of velocity and 
physical extension are shown in Figures 8 and 9. ARMF-MM is able to maintain a low level of 
velocity RMSE in the entire process and performs better than the others in physical extension estimation. 
The physical extension RMSE of ARMF-MM is significantly reduced when the target maneuvers. 
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Figure 8. RMSE of velocity in simulation scenario 2. 
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Figure 9. RMSE of physical extension in simulation scenario 
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6. Conclusions 

Bayesian EOT using a random matrix is a very effective and simple approach to incorporate 
physical extension in the OT framework. The Bayesian RMF approach was reviewed first in this 
paper, and certain improvements were made. Afterward, model parameter adaptive approaches were 
derived for both extension dynamic evolution and measurement noise based on the analysis of SPD 
matrix decomposition. The proposed adaptive RMF can be integrated into an IMM framework to 
further improve estimation performance. The validity of ARMF and ARMF-MM was verified by the 
simulation results. ARMF is able to change model parameters adaptively when the extended object 
maneuvers, thereby resulting in low physical extension estimation error. 
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